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On a basis of a two-length scale description of hydrophobic interactions we de- 
velop a continuous self-consistent theory of solute-water interactions which allows 
to determine a hydrophobic layer of a solute molecules of any geometry with explicit 
account of solvent structure described by its correlation function. We compute the 
mean solvent density profile ra(r) surrounding the spherical solute molecule as well 
as its solvation free energy AG(r). We compare the two-length scale theory to the 
numerical data of Monte-Carlo simulations found in the literature and discuss the 
possibility of a self-consistent adjustment of the free parameters of the theory. In the 
frameworks of the discussed approach we compute also the solvation free energies of 
alkane molecules and the free energy of interaction AG(D) of two spheres of radius 
R separated by the distance D. We describe the general setting of the self-consistent 
account of electrostatic interactions in the frameworks of the model where the wa- 
ter is considered not as a continuous media, but as a gas of dipoles. We analyze 
the limiting cases where the proposed theory coincides with the electrostatics of a 
continuous media. 

I. INTRODUCTION 

We consider the current work as a step towards the development of a continuous thermo- 
dynamic theory of solute-solvent interactions. The construction of such a kind of a theory 
is of extreme importance for biochemical applications connected with the reliable and fast 
determination of the protein-ligand binding constants in solvent (water) with the precision 
comparable to the accuracy of explicit molecular-dynamical simulations. 

There are few kinds of approaches which take into account the effect of fluctuating media 
on interactions between solvated molecules. Tentatively we can split them into the following 
three groups: 

i) "Explicit approaches". Generally such methods are based on numerical simulations 
taking into account motion of many water molecules. Such approaches are very time con- 
suming because they demand huge computational resources. 

ii) "Implicit approaches". These: a) either exploit the heuristic concepts of the "surface 
area accessible by the solvent" and of the "width of the hydrophobic layer" , or b) make some 
effective renormalization of "bare" (i.e. vacuum) parameters of interacting molecules: in the 
frameworks of these approaches the effect of water on thermodynamic properties of dissolved 
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substances is computed on the basis of some additional data on solubility, evaporation etc— 
see, for example, (H 0,111 • 

iii) "Intermediate approaches". The works belonging to that group (including our one), 
offer attempts which, on one hand, get rid of explicit numerical simulation of water molecules, 
and, on the other hand, preserve an effective account for water. The development of these 
methods on the basis of a clear understanding of a physical origin of solute-water interactions 
could, at least, give a better control for the parameters used in implicit approaches. 

Most commonly, in the "implicit approaches" it is assumed 0, H| that the free energy 
of solvation is proportional to the solvent accessible surface area of a solute molecule. It 
is obviously true for solutes much larger than a water molecule. However, for the solutes 
of sizes comparable to the size of water molecules, as well as for the solutes of complex 
geometry, the question of the determination of the corresponding hydrophobic layer, being 
of great importance, is often resolved by means of purely empiric conjectures. For example, 
it is supposed usually that the hydrophobic layer is equal to the surface covered by a center 
of a water molecule "rolling" around the solute surface. Some modifications of the method 
deal with the so-called "molecular surface area" - see, for example, [jj. Another method of 
implicit account for the solvent in the so-called " Gaussian approximation" has been proposed 
by 0. 

Among the theoretical attempts to develop a constructive theory of solvent-solute inter- 
actions, the following groups of works (to our point of view) deserve special attention: 

1. The detailed computation of the solvent density-density correlation function in the 
presence of dissolved substances (with or without electrostatics): a) either introducing 
weighted densities and explicitly designing the form of the density functional which is en- 
tirely dictated by physical constraints or b) by the construction of an appropriate "bridge 
functional" (i.e. closure relation) for the correlation function satisfying the hierarchy of inte- 
gral equations 0. On the basis of the correlation function all the thermodynamic quantities 
can be easily determined. 

2. The computation of the averaged profile of the solvent density described by the fluctu- 
ating field in presence of solute molecules (with or without electrostatic interactions). The 
solvent structure on small scales is taken into account by its bulk correlation function, while 
the effect of solute solvation is considered by forcing the total solvent density to be zero in- 
side the solute molecule (see, for example, [lOj] and references therein). Minimization of the 
corresponding free energy allows one to compute the density and the solvation free energy. 

3. The "Scaled Particle Theory^SPT) Q [lj and the "Information Theory" (IT), 



developed in the series of papers [13]. Both these approaches, SPT and IT, exploit the 
relation between the probability to find a cavity of a given shape containing some number of 
particles, with the solvation free energy of a solute molecule of that shape. The probability 
of a cavity formation in an ensemble of fluctuating particles (water molecules) is computed 
on the basis of semi-empirical combinatorial arguments. 

Our solvation model borrows its ideas in the second group from the list above, namely 
from the series of works 0, 0, 0] where a two-length scale description of hydrophobic 
interactions was developed. In this approach the solvent density is decomposed in two 
components, the slowly varying field describing the mean solvent (water) density, and the 
short-ranged density fluctuations. We believe that such an approach is optimal from different 
points of view: on one hand it is physically clear, being "semi-microscopic", and on other 
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hand it can be used as a constructive computational tool of account for water, much faster 
than the corresponding explicit approaches but without the loss of the precision. So, our 
aim is to construct a theory of solute-water interactions which would allow to determine 
the hydrophobic layer of solute molecules of any geometry with account of water structure 
described by its correlation function. The parameters of such a theory will be determined 
by testing them for solutes with simple geometries and comparing these cases with the 
experimental data. Briefly, the main ingredients of our work are as follows: 

i) We consider water as a continuous inhomogeneous media and describe it by a fluctuating 
continuous density field; the discreteness of the water structure on small scales is taken into 
account through the fluctuations controlled by the water correlation function x i n the bulk; 

ii) The solvation free energy of the substance situated in water is obtained in the self- 
consistent mean-field approximation based on slightly modified Ginzburg-Landau-Chandler 
functional 3 HI ■ 



Our work is devoted to the calculation of the solvation free energy and of the free energy 
of interactions of solute molecules of any arbitrary shapes. 

Specifically, we start with a two-scale Ginzburg-Landau Hamiltonian and perform an 
analytical integration over the short-ranged fluctuations. This provides a single-variable 
Hamiltonian yielding an integro-differential equation. The solution of this equation ensures 
the equilibrium mean density profile of the water, and the corresponding value of the Hamil- 
tonian gives the solvation free energy. We pay a special attention to the physical clarity 
and consistency of all steps of derivation of our model. We discuss the hidden dangers on 
the way and compare our approach to the Monte-Carlo numerical simulations carried out 
in 0| • We determine in a self-consistent way the width and the geometry of a hydrophobic 
layer surrounding the solute molecules. So, we believe that a very important auxiliary goal 
of our work concerns the prediction of the relevant number of water shells (the "solvent 
accessible surface area") around a solute molecule. 

Summarizing the said above, the advantages of the fluctuational approach with respect 
to other phenomenological theories of " implicit water" are as follows: 

a) There is no difference in computations of the solvation free energy and of the free 
energy of interactions of solute molecules; 

b) The detailed structure of the solvent is taken into account on small scales; 

c) We do not need any special determination of the "hydrophobic layer" by "rolling" the 
water molecule around the solute surface. The width and the geometry of a hydrophobic 
layer is computed on the basis of: a) the true geometry of a solute molecule, and b) the 
supposition about the solute-selvent and solvent-solvent interactions. 

d) The theory contains small number of adjustable parameters, which in turn can be 
determined from the independent experiments (either numerical, or real). 

Our paper is organized as follows. In the Section II we describe the model and derive the 
basic equations for the equilibrium and the full densities, and for the solvation free energy, as 
well as we solve numerically the obtained equations, compare them to the numerical Monte- 
Carlo simulations of the work and propose the way of tuning the adjustable parameters 
of the Hamiltonian; in the Secion III we compute the solvation free energy of alkanes, as well 
as the free energy of interaction of two spherical solute molecules; the results of the work 
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are summarized in the Section IV, where we also derive the equations for electrostatics in 
fluctuating dipolar environment. 

II. THE MODEL 

In this Section we present the self-consistent description of solvation which takes into 
account simultaneously: i) the smoothly varying average density field n(r), and ii) the water 
structure on small scales described by the density-density correlation function x(r,r'). 

Following the general scheme of the works 16], we describe the solvent by two fields: n(r) 
- the average density varying on long distances, and u;(r) - the small-scale fluctuating field 
controlled by the correlation function of the water x(r,r') = (u;(r)c<;(r')) in the bulk. The 
solvent cannot penetrate into the volume occupied by the solute molecule hence the total 
solvent density p(r) = n(r)+u(r) is nullified inside the solute: p(r) = for all r G v- in , where 
f in is the volume occupied by the solute molecule. The fact that we force the total solvent 
density p(r) = n(r) + ui{r) to be zero inside the solute, results in an effective interactions 
between n(r) and uj(t). In general, one can permit also the direct coupling between n(r) 
and u;(r) everywhere in the solute. 

The thermodynamic properties of the solvent (water) in the bulk are obtained using the 
self-consistent Ginzburg-Landau (GL) free energy functional which describes the liquid- 
vapor phase transition. The average solvent density n(r) plays the role of an order parameter. 
The interactions between long- and small-scale fields n(r) and o>(r) can enforce the "drying" , 
shifting the liquid-vapor transition. Minimizing the corresponding free energy functional, 
we compute the desired density profile and the solvation free energy. 

It should be noted that the fully microscopic description of the solute-solvent interactions 
certainly does not require introducing of any scale separation — everything is described by the 
single microscopic field p(r). However in such a description the interactions between water 
molecules should also be described on the microscopic level. Such kind of theory has been 
developed in the works 0]. It seems to us that due to the heavy machinery, the application 
of such approach to the solvation and interactions of objects of complex architecture in water 
is yet rather restrictive for practical purposes. So, the nature of two different scales in the 
approach chosen in our work, is the consequence of the mean-field (GL) description of the 
water-water interactions. 

Before going into details let us point out what is different in our consideration with 
respect to the original two-length scale theory of Chandler et al. [1(1 Eil [Til fl6j] . 

i) We do not divide the space in "cells" for the large-scale field n(r) (as in [lol|). We 
describe the system by two fully continuous densities. The equilibrium density is determined 
selfconsistently without additional averaging over small-scale fluctuations. We believe that 
this would result in better description of solvation of objects of rough shape and of the free 
energy of interactions between solutes located close to each other. 

ii) We can tune the potential W(n(r)). Currently W(n(r)) describes the behavior of the 
system exactly at the vapor-liquid phase transition, but our consideration can be straight- 
forwardly generalized to other situations. 

iii) We can consider the charged objects in the frameworks of the "discrete electrostatics", 
where the density n(r) becomes the density of water dipoles. The set of basic equations for 



5 



the density field /;(r) and the electric field 0(r) is obtained in Section III, while in the 
forthcoming paper [17[ we numerically solve the derived system of equations. 

iv) We can incorporate in the theory the direct Van-der-Waals solute-solvent potentials 
and consider the smooth solute-solvent boundary (see the Section ITTT1 for more details). 



A. The Hamiltonian 



We describe the solvent in absence of any solute molecules by the following two-length 
scale Hamiltonian 

H [u>,n) = ^J W (r)f 1 (r,r>(r')tt' + J {- (Vn(r)) 2 + W(n(r))} dv (1) 

where n(r) is the smoothly changing (average) solvent density (also considered as the "order 
parameter" of the theory); u(r) is the field corresponding to the short-ranged density fluc- 
tuations; x(|r — r'|) is the solvent correlation function in the bulk; a is the phenomenological 
parameter which requires the microscopic determination — see the discussions below; and the 
self-consistent potential W(n(r)) is chosen in the common form of the standard Ginzburg- 
Landau (GL) expansion for the order parameter n(r) as the fourth-order polynomial: 

b 

W{n) = -(n-n^ 2 (n - n 2 ) 2 (0 < n x < n 2 < 1) (2) 

where n\ and n 2 are the values of the order parameter n in the vapor and liquid phases 
correspondingly (below we set, if not specified, n\ = 0) and b is the coupling constant which 
in combination with the parameter a defines the surface tension (see fjlOj) below). 

The Hamiltonian consists of two (still decoupled) parts. The first term describes the 
non-local small-scale Gaussian fluctuations of the solvent. These fluctuations extend to the 
distances controlled by the solvent density correlation function in the bulk %(|r — r'|), 

X (|r-r'|) = (c(rV(r')) (3) 

The second term describes the fluctuations of the self-consistent potential (0) in the system 
which could exhibit a phase transition. In our case this is a liquid-vapor (i.e. "drying") 
phase transition. 

The choice of a functional form of the potential (0) is rather arbitrary being to some 
extent a question of a taste. However it should obey some basic requirements which are 
imposed mainly by physical reasons. As it has been mentioned already, the Hamiltonian 
(JTJ) with the potential Q should describe the vapor-liquid (i.e. "drying") phase transitions. 
This should be valid only in the vicinity of the solvated object, thus manifesting the nature 
of a hydrophobic effect. The presence of a drying transition near the solute surface shifts 
effectively the parameters of the bulk media to the region of the liquid-vapor coexistence 0] . 
The general form of the Hamiltonian describing the liquid-water coexistence is well known 
- it is the two-well potential usually parameterized by the GL expansion — see Fig[T] The 
parameterization being rather arbitrary, might be taken not only in the polynomial form as 
in the GL theory, but also in any form that correctly describes the region of phase coexistence 
under conditions close to normal. We believe that the GL expansion is well enough to reflect 
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the major features of the physical picture. The choice of a more sophisticated form of the 
potential (J2J) entering in the large-scale part of the Hamiltonian (Q) is the possible way of 
the refinement of the model in future within the frameworks of the mean-field approach. 

In FigQ]for generality we show broader class of potentials Win) than described by Eq.(J2J), 
allowing the nonsymmetric shapes (see FigO^)- The profiles in Fig^i,b are drawn for n\ = 0. 
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(a) (b) 

FIG. 1: Sample plots of various shapes of the potential W{n). The case (a) corresponds to and 
is analyzed in this paper. 



The potential W(n) in the GL expansion depends upon four parameters (a, b, n 1; n 2 ), 
two of them (the densities n\ and n 2 ) are fixed in each phase. The parameters a and b are 
considered as the free adjustable coefficients of the theory. Supposing slightly more general 
form of the potential 121 Win) = Win) + u(n — ni), shown in Fig^o, we may consider 
the relative difference (RD) between the energy densities in each phase, i.e. the difference 
between the values at the minima of W{n) (see FigO^) as an extra adjustable parameter. 
It is easy to see that this difference is closely connected to the asymptotic behavior of the 
solvation energy at large solute sizes. Actually, the RD contributes to the volume-dependent 
part of the solvation energy. It is obvious that the term which scales linearly with the volume, 
tends to |PV^ iutc as the solute volume Volute increases (here P is the liquid pressure). 
However the simple quantitative estimation shows that on the nanometer scale the volume 
contribution is negligible in comparison with the surface contribution to the solvation energy 
|lo| . Thus we may neglect the difference in minima of the potential shown in FigO^ and 
consider everywhere further the potential symmetric upon the phase change — see Fig^i. 
Just such functional form of the potential is proposed above in (J2J). 

We can find the equilibrium densities minimizing the Hamiltonian (0) with respect to the 
(still decoupled) fields n(r) and p(r): 

J x(r - r>(r') dr' = (4a) 

. . , 5W(n(r)) 
-oAn r + f V, = 4b 
on(r) 

In absence of any interactions between the fields u and n, the solution of (|4aj) postulates the 
homogeneity of the liquid uj{r) = const, while Eq. (j4b|) describes: a) the plain vapor-liquid 
interface supplied with the boundary conditions n(— oo) = n\ = (n(oo) = n%) in the one- 
dimensional case, and b) the density profile near the surface of extended macroscopic object 
in the two- and three-dimensional case with boundary conditions n = at the solute surface 
and Vn = for |r| — ► oo. Everywhere further in this Section for simplicity we consider the 
spherical objects only (the solvation of objects of other shapes is discussed at lengths of the 
Section III). 



Qualitatively these solutions can be analyzed by computing the first integral of motion 
(as in the classical mechanics [lij). Writing the first term in Eq.(|4b|) in 3D as An(r) = 
n"(r) + -n'(r), multiplying ()4b|) by n'{r) = and integrating it from r in (the radius of 
the sphere) to oo, we obtain: 
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-(n'(r)) 2 dr + 
r 



dW(r) 
dr 



dr = 0; 



n"(oo) = (5) 



Hence, we get from (|SJ) the following relation: 



(n'(r))' 



n{r-)ydr + -[^(n(r in )) - W(n(oo))} 



(6) 



Having the equation © and the exact functional form of the potential (J2J) in hand, we know 
the value of n'(r) at the point r = r in which determines the qualitative behavior of the 
density profile. In our case the second integrand vanishes, hence the density monotonically 
increases from n\ at the point r in to n 2 for |r| = oo. In practice, the liquid density reaches 
quickly the bulk value at some distance outside the solute surface. The analysis of an 
asymmetric potentials with RD can be performed within the same formalism. 

The sense of the parameters a and b becomes transparent in the limit of large spherical 
solute molecules r- m ^> 5, where 5 is the characteristic size of the liquid-vapor interface (see 
below). Substituting (j2J) into (}4b"|) and replacing the interaction of the field n(r) with the 
solute molecule by the boundary condition, we get 



-dAn(r) + 2b(n(r) — n avT )(n(r) — ni)(n(r) — n 2 ) = 
n(r = r in ) = n t 



(7) 



where n R 



n 2 + ni 



Supposing in ((7|) that r ^> 5, we can expand the 3D Laplacian 

d 2 
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d 2 2d 
dr 2 r dr 



r><5 



dr'" 
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and obtain the effective ID equation 

—an"(r) + 2b (n(r) — n avr )(n(r) — ni)(n(r) 

^('"in) ^"avr 

whose solution reads 
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n(r\ 



n 2 - nx 



tanh 
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n 2 + ni 
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n 2 - rix 



(8) 



(9) 



It is clear that the symmetric potential of GL form leads to the one-parameter theory. This 
parameter 5 fixes the width of the phase separation interface. The solvation free energy G 
per the spherical solute area 47T7f n in the limit of large molecules can now be expressed as 
follows 



G 
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where the coefficient 

Vab 3 

o = —^-{n 2 - ni) 
o 

has a sense of a surface tension. Note that for < — <C 1 in the leading order of the 

^*in 

expansion (fTUj) the true boundary condition n(r in — 5) = cannot be distinguished from the 
approximate one, n(r in ) = n avr . 



B. The correlation function 



The important component of our model is the bulk correlation function, x(|r — r'|), of 
a solvent. It appears as an input into the theory. Recall some basic facts concerning the 
function %(|r — r'|). We have by definition: 

X(|r - r'|) = nS(r - r') + n 2 h(\r - r'|) (11) 

where n is the mean solvent density in the bulk, h(r — r') = g(r — r') — 1 and the function 
g{ I r — r' | ) has very clear physical meaning: 



n 2 



rMr) = y- 1 <^5(r-r l + r,)> (12) 



The function g(r) may be determined from experimental data, or numerical simulations, or, 
say, from some selfconsistent integral equation with corresponding bridge functional |l9j . 

The following relation is the direct consequence of the definition of x{ r ~ r ') : 

n J (g(\r - r'|) - l) dr' = Ann J h{r)r 2 dr = -1 (13) 
The inverse correlation function X\2 ( r ' r ') i s defined as follows 

f dr" X ^(r, r") X (\r" - r'\) = 5(r> - r) for 



r G v in 

Xm( r . r ") =° for r G v out 



(14) 



Let us pay attention that the integration in ()14j) is spread only to the volume v- m occupied 
by the solute and hence the function (r, r') is not translationally invariant. 

We have used in our work two different correlation functions: i) the correlation function 
of water obtained from the numerical molecular-dynamical simulations 2ol |. and ii) the 



correlation function of hard spheres constructed on the basis of solution of self-consistent 
integral equation completed by the Percus-Yevic bridge functional ji^, • The functions 
h(\r— r'|) = g(|r— r'|)— 1 of water and of hard sphere liquid are shown in FigJ21 Both functions 
correspond to the dimensionless bulk density n 2 = 0.7. Some approximate equations of the 
theory of liquids in the statistical thermodynamics of classical liquid systems one can find, 
for instance, in the review article [l9| . 

Working simultaneously with two different correlation functions (of water and of hard 
spheres) we able able to check the sensitivity of the physical observable quantities (as the 
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FIG. 2: Correlation function of water vs Percus-Yevic correlation function of hard spheres. The 
distance is measured in dimensionless units r/£ (£ = 2.78 A), the dimensionless bulk density is 
n 2 = 0.7. 

density and the free energy) to the details of the correlation function, as well as to verify 
the convergence of our numerical scheme: far away from the solute surface the results of 
computations should not be sensitive to the particular choice of the correlation function. 

The correlation function defines the natural scale, £, in the theory. In what follows we 
fix £ = 2.78 A as the distance till the first peak in the correlation function of water. 




C. The partition function 



Now we turn to the solute-solvent interactions. To take into account the effect of the 
solute solvation, we demand the total density p to be zero inside the volume v- m of the solute: 

p(r) = n(r) + w(r) = (15) 

where r G v- in . The direct interaction between the short- and the long-range density fluctu- 
ations in the solvent can be written in the simplest possible way 

V[u,n} — c J lu(t) n(r) dr (16) 

where c is the coupling constant. In what follows all the analytic derivations are performed 
for generic c, however the numerical solutions of the corresponding equations are presented 
for c = only. The case c^Owe leave for the future consideration. 

The solute-solvent partition function is obtained as follows 

Z= f V{u{v)} J] 5[p(r)]e- n (17) 

J re V in 
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where 

H = n [u,n] + V[uj,n} + J drtp{r)p{r) (18a) 

p(r) = n(r)+w(r) (18b) 

Let us rewrite the constraints imposed by the Dirac 5-function in (fTTj) using the functional 
Fourier transform 

II 6 ^ = J2^W I WW} ex p{ z / PtotKr)dr| ( 19 ) 
where Af is the number of points in the product f] (...), and we have denoted D{ip(r)} = 

re Win 

Yl dip{r). Now we can represent the partition function Z[y>(r)] in (fTTjl as a functional 

r£v in 

integral over an auxiliary field ip(r): 

Z = '~- 1 I V{uj(r)}V{iP(r)}e- s ^ n ^ 

f (20) 
S[u>, n, tp] — H[ui, n] + i / p(r) vp(r) dr 



where H is some physically irrelevant constant which defines the normalization of the parti- 
tion function Z. 

The partition function Z with the action S[u,n,ip] defined in (J2U|) is the key object for 
our computations. Such quantities as the total, p(r), and the averaged, n(r), densities, the 
solvation free energy, AG, the free energy of interactions of different solute molecules in the 
solvent, etc. can be calculated on the basis of the generating function Z. 



D. The mean density 

To define the mean (large scale) density n(r), let us require 

SS[n(r)} 



8n(r) 



(21) 



In (121 J) all other fields except n(r) are supposed to be already integrated over. This leads 
to the functional equation which can be written explicitly as a set of two coupled integro- 
differential equations relating the inner (r G v in ) and outer (r G v out ) regions of the solute 
molecule: 

. , , 5W(n(r)) , . 
-oAn r + ; \ \ ~ 2cn(r)- 
on(r) 

c f dr' I dr" X r\ r , r')x(|r' - r"|)n(r") + / dr' X r n l {r, r')n(v') = for r G v in 

(22a) 
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— c 2 / dr'x(|r — r'|)n(r')+ 

•'"out 



and 

. , > 5W(n(r)) 

c 2 [ dr' I dv" [ dr"' x (\r-r'\) x r\r',r") x (\r"-r"'\)n(r"')- (22b) 

Jv in Jv in Jv out 

c dr' r"x(\r - r'\)x^(r', r")n(r") = for r e v ont 

Jv in Jv in 

In the current paper we pay attention to the case c = only, i.e. when the direct coupling 
between the short- and long-ranged fluctuation is absent. In that case the profiles of the 
mean density n(r) for spherical solute molecules of different sizes are given by the solutions 
of the simplified integro-differential equation 



- aAn(r) + 6W ^^ + f dr' X ^{v, r')n(r') = for 
on[r) J 



r G v- m (23a) 



-aAn(r) + 5W( >^ =0 for r e v out (23b) 
on(r) 

We solve (|23a|) - ()23b|) numerically by using the simple technical trick which allows to 
increase essentially the speed of the computations. Namely, we have observed that the most 
CPU time is spent for the calculation of x^n( r ^ v ') on the basis of (JHJ). To get rid of this 
part of computations we can multiply the ()23a|) by x(\ r ~ r '\) an d integrate over the whole 
space v to t (v to t = ^in U f ou t), with the subsequent application of (ftl| . Finally we arrive at 
the following set of equations: 

n(r) + [ dr' X (\r - r'|) (-aAn(r') + 5 J^}^J1 ] = for r e v in 

. , SW(n(r)) 
-aAn(r) + \\ =0 for r e w out 
on r 



The generic form of the mean density profile n(r) for the spherical solute molecule is 
shown in FigEJfor different (still arbitrary) values of the constants a and b of Hamiltonian. 
The densities of the vapor (rii) and bulk (722) phases are as follows: n\ = 0; ri2 = 0.7. 
The numerical solutions of (J24j) are obtained with the correlation functions of hard spheres 
(dashed line) and of water (solid line). The typical plots are shown in FigJ21 

The presence of the regions with the negative values of the mean densities in the vapor 
phase should not confuse — remember that only the full density p(r) is a conserved value and 
should be zero inside the solute molecule, while the average density n(r) is determined self- 
consistently. The existence of the regions with negative values of the long-ranged averaged 
field n(r) signifies that in these regions the fluctuations of the short-ranged field p(r) are 
not symmetric, such that locally (u) 7^ 0. Such regions exist only inside the solute volume 
f in and in the "physical" volume t> out (i.e. outside the solute) the mean density is always 
positive. 
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FIG. 3: Typical profiles of the mean density n(r) around spherical solute molecules of different 
radii Ri n . 



The full density p(r) = n(r) + uj(r) which minimizes the functional (J2(J|) satisfies the 
following set of equation 

p(r) = 0, for r G v- m 

p(T)=n(T)-f dr' I dr" X (\r-r'\) X r 1 (r',r")n(r")+ 

c dr' dr" dr"'x(\r -r'\)xrn(r',r")x(\r -^"Dnir'^- 

JVin Jv in Jv in 



dr'x(\r — r'|)?2(r') for r G f out 

(25) 

The equations (|25|) with c = again can be rewritten in the form which excludes the 
function ^^(r,! - ') from the corresponding expressions. 



E. The solvation free energy 



The free energy of solvation AG is defined as the energy necessary to transport the 
solute molecule from its environment in the solvent to the vacuum. The partition functions 
of solvent samples of sufficiently large volume Solvent containing the solute molecule inside 
can be straightforwardly written on the basis of (fTTj) : 



J solvent — solute 



/ 2>{"(r)} II W r )l e ~ W ( 26 ) 

•^solvent rei> so iutc 

where 7i is given by (J18a|) . Conventionally the solvation free energy AG is written as follows 

a i ^solvent— solute / 

AG = -ln — (27) 

•^solvent only 
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where Solvent only is the partition function of the pure solvent (without the solute molecules) . 
At first glance it seems naturally to write Solvent on i y simply as 



J solvent only 



V{co(r)}e- n (28) 

^solvent 



However the expression (|28|) being used in (|27jl leads to the divergence of the solvation free 
energy AG in practical computations based on the approach developed in ([T7 ]) -pi j) . The 
formal reason for such a divergence deals with the occurrence of uncompensated infinite 
product of Gaussian integrals (|20*j) in the ratio Solvent-solute /Solvent only in (EH)- The physical 
origin of such a divergence is due to the mixture of different statistical ensembles associated 
with the partition functions (j2*U|) and (|25jl. Namely, when writing Solvent-solute as in (|2*U|) 
and imposing the ^-function constraint, we fix some particular value of the field p within the 
solute volume v- m what means that (|2fi|) is the partition function of the canonical ensemble 
(with respect to the density p). At the same time, for Solvent only written in the form (J25J) . 
we allow any density of the field p inside the solute volume. Hence, (J2H|) is the partition 
function of the grand canonical ensemble. 

The regularization of eq. ()28|) is based on the probabilistic consideration of the solvation 
free energy Rewrite AG in (j2*7|) as follows 

^--sS&o (29) 

where Z V (N) is a partition function of N molecules inside the volume v- m . The continu- 
ous analog of the last expression appears when passing from N to the average number of 
molecules n = N/v- m . Simultaneously we require the net density p to be equal ft inside 
f; n and replace the summation over N by the integration over n. Now we can rewrite the 
normalization partition function Solvent only in the following "regularized" form 

/oo p 
dh I V{u(r)} Yl 5 [p(r) - n] e~ n (30) 
•^solvent rev solutc 

Instead of this form we use an approximation 



J solvent only 



/ V{u(r)} Hp(r)-n]e- n (31) 

" ^solvent — 



resolute 



obtained from (j3*Uj) by the stationary phase integral calculation technic. Thus, the definition 
(I27jl with the partition functions given by (j2Hj) and (|3*T|) is very natural, justified physically 
and does not contain any divergences. It is important to notice that the computation of 
the equilibrium density based on Solvent only given by (JBTij) is the bulk density, for which the 
value of corresponding effective Hamiltonian is strictly zero. 

The solvation free energy AG for any solute molecule is computed on the basis of the 
equilibrium density profile n eq (r) = n(r) calculated from Eqs. (|22a|) - (|22b|) . Substituting 
n e q(r) back into the Hamiltonian (|18aj) and taking into account (|14|) . we can rewrite the 
solvation free energy G = H.{n eq } as a sum of three terms: 

AGgolv^in) = AG C orr + AG sur f + AGhit (32) 
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where 



A /"> 

AG sur f 
AGi nt 



5 / ^cq(r)x in 1 (r,r / )n cq (r / ) drdr' 



(Vn eq (r)) 2 rfr 



W (n eq (r))dr 



(33) 



As in the computations of the mean density n(r), here we consider the case c = only. 
In that case ()32j) can be rewritten in the following form 



AG so iJv in ) 



jxr / \ n8W\ 
W(n) - I dr 



2 5n 

The contributions to the normalized solvation free energy 

AC7 so l v 



(34) 



AG, 



solv 



Arcrt 



(35) 



of a spherical solute molecule of radius R are shown in Fig|l]for each term {AG COTT , AG sur f, 
AGint} = {AG C orr, AG sur f, AGi nt }/ (47i\fif n ) separately, as well as for the sum AG so i v = 
AG C orr + AG sur f + AGi nt . The diameter of the solvent molecule in FigE|is set to 1. 





FIG. 4: Typical behavior of the solvation free energy for different choice of parameters a and b 
in the Hamiltonian (left); Contributions to the normalized solvation free energy AG so i v (i?; n ) of a 
spherical solute from AG corr , AG sur f and AG; n t for a = 9, b = 30 (right). 



The solvation free energy of a spherical solute molecule in absence of electrostatic in- 
teractions [U, EH, is proportional (as expected) to the volume of solute molecule for 
sufficiently small sizes i?i n (of order of the correlation length in the solvent) and tends to be 
proportional to the surface area for large R{ n . We clearly see the non-monotonic behavior 
of the solvation energy upon the size of a solute molecule j3^| reported in some works (see, 
for example, [lol]). 
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F. Adjustment of the parameters of the Hamiltonian 



The free parameters a and b of the Hamiltonian (JIJ-(J2|) are chosen from two basic re- 
quirements: 

i) The theoretically computed normalized solvation free energy AG solv = AG solv / (4tcR 2 ) 
of a spherical solute molecule of a radius R (see (j35| ) reproduces the corresponding depen- 
dence AG so i v (R) obtained in the Monte-Carlo simulations of [Toj | : 

ii) The theoretically computed full density p(r\R) for few sizes R of spherical solute 
molecules (|25jl reproduces the behavior g(r) of the correlation function in vicinity of the 
spherical solute of radius R extracted from the Monte-Carlo simulations [nil ]. 

The corresponding results are shown in FigJ3]for the following choice of the parameters: 
a = 9; b = 30. 



1.2 




.0,1 I I I I 

0123456789 10 11 0123456789 10 11 12 



r[A] r[A] 



FIG. 5: Solvation free energy (|34[) of a spherical solute vs Monte-Carlo simulations [lfj (left); Pro- 
files of the full density p(r) around spherical solutes of different radii R vs Monte-Carlo simulations 
flol | for the correlation function g(r). 



The dashed and solid lines represent the results obtained with the correlation functions of 
hard spheres and of the water correspondingly. The solvation energy AG so i v is represented in 
FigEJin dimensionless units AG so i v /7, where 7 is the surface tension of the flat vapor-water 
interface. The numerical value 7 ~ 72 mJ/m 2 is taken from (lp| . 

The solvation energies AG so i v and AG so i v computed directly using (J3~4"j) and (|33j) can be 
converted to the proper units by means of the following scaling coefficients: 

AG so iv, [dimensionless units] x 0.2479 = AG so i v , [kcal/mol] 

0.2479 ~ (36) 

AGgoiv, [dimensionless units] x — = AG so i v , [kcal/(mol A 2 )] 

(2.78) 2 



These transfer coefficients remain unchanged in all further computations. In Fig0]and FigEl 
the radius of a solute molecule is measured in Angstroms. 
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III. RESULTS AND DISCUSSIONS: SOLVATION OF OBJECTS OF VARIOUS 

SHAPES 



The current mean-field-type theory of solvation is applied to computation of solvation of 
neutral molecules (alkanes), as well as to the computation of the free energy of interactions 
of separated objects (spheres). 



A. Solvation of alkanes 



The parameters of the theory and the scaling coefficient are adjusted to satisfy the Monte- 
Carlo simulations of solvation of hard spheres and the corresponding profile of the solvent 
density near the solvated object (see the previous Section) and are not tuned anymore. In 
particular, we use the same parameters to compute analytically the solvation free energy of 
alkanes. As it is shown below, we find very good coincidence of our computed values with 
the predictions of the " Scaled Particle Theory" which uses the parameters tuned especially 
to alkanes ul 



The mentioned coincidence needs some elucidations. There is a viewpoint [27|, |28|, |2£ 
supported by numerical computations that one can split the interaction between the solute 
and the solvent into two parts: i) the free energy, AG cav , of a "cavity formation", and ii) the 
dispersion (attractive) part of the Van-der-Waals interactions, AGdisp- The numerical meth- 
ods involving the "thermodynamic integrations" allow to compute both the contributions, 
AG cav and AGdi sp separately These contributions can be also computed separately in 
the frameworks of the "Scaled Particle Theory" (SPT) developed in The correspond- 



ing description of alkanes has been undertaken in [25] where the authors have reported very 
good quantitative agreement of the sum AG cav + AGdisp with the experimentally measured 
solvation free energy. 

The model described in the previous sections of our work considers the solute as an ob- 
ject bounded by hard walls and hence takes into account only the "cavitation" part of the 
solvation free energy. The attractive part of the Van-der-Waals interactions is not yet taken 
into account. However we do not see any principal obstacles in adding the dispersion part 
of solute-solvent interactions to our model. Moreover, we can "smear" the (^-functional 
constraint in (|19|). releasing more realistic form of the Van-der-Waals potential. The cor- 
responding computations are in progress and will be reported in a forthcoming publication 

u 

The comparison of our predictions for the cavitation part of the solvation free energy of 
alkanes with the corresponding contribution AG cav extracted from the papers [25| and [26| 
is shown in the table below. 

Let us remind once more that in our model the parameters of the Hamiltonian a = 9, 
b = 30 and the scaling factors (jHUj) are tuned on the basis of the Monte-Carlo simulations 
of solvation of hard spheres without any additional adjustment to alkanes. The volume v- in 
occupied by the alkane molecules we compute using the following prescription. We take 
the coordinates of centers of all atoms from a standard database, then we draw the spheres 
with the effective radii 3.5 A for all carbon atoms (the length of the C-0 bond in water) 
and 3.05 A for all hydrogen bonds (the length of the O-H bond in water) and find v- m as 
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1 K^TIPK 


/\(-r Kfal/iTiol ( our monpli 


A (7 kcal/mol (from [2511 


AG kcal/mol (from [2611 




o.oy 


O.Dl 


O.OD 


P IT 

^2*16 


( .04 


< .00 


1 . 10 


C3H8 


9.07 


9.13 




C4H10 


10.52 


10.8 




C5H12 


11.96 


12.8 




C6H14 


13.24 


14.8 


14.22 



TABLE I: Comparison of the "cavity formation" contribution, AG cav , to the solvation free energy 
of alkanes computed on the basis of our model with the corresponding data extracted from the 
papers and 0. 



the union of the corresponding spherical volumes. The corresponding numerical values of 
C-0 and O-H bonds are extracted from the Monte-Carlo (MC) simulations of the sytemem 
consisting of 216 water molecules and one methane molecule and subsequent computations 
of the pairwise correlation functions [C(methane)-0(water)] and [H(methane]__p(water)], 
*30l |. The same parameters can be extracted using the force field developed in |3l| . 



As it has already been mentioned above, the dispersion part of the Van-der-Waals (VdW) 
interaction can be easily taken into account. In the presence of the attractive contribution 
to the VdW interactions, we should add an extra term to the mean-field potential (JTJ). 
where 

f^disp = J U(r) n(r)dr 

and U(r) is the attractive part of the VdW interactions. 

In the presence of the potential U(r) the equations ()23a|) - (|23bJ) for the equilibrium mean 
density are modified in the following simple way: ()23a|) reminds without changes, while the 
equation (J23bj) acquires an additional contribution from the attractive part of the potential 
U(r) in the region r G v out : 

_ aAn(r) + f2^M> + u(t) = (37) 
on{r) 



B. Solvation of cylinders and the free energy of interactions of two separated 

spheres 



In this section we consider the cavitation part of the solvation free energy of molecules 
of cylindrical geometry, as well as of the free energy of interaction of two hard spheres in 
the fluctuating media (i.e. in water). Schematically the studied systems are shown in FigEJ 
The purpose of this section is two-fold. 

1. First of all, we compute the depth of the "skin layer" of the solvated object. The naive 
computation of the solvation free energy (as well as of the free energy of interactions) of large 
molecules like ligands and proteins demands essential computational resources. Actually, 
according to ({21} one has to take into account the interactions of the fluctuating density 
field n(r) with all points inside the solvated object v in . However it is clear that due to the 
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(a) 



D 

(b) 



FIG. 6: Schematic pictures: of a solvated cylinder ended by two semi-spheres (a), and of two 
spherical solute molecules (b). 

hard wall constraint the fluctuating field n(r) cannot penetrate deep inside the body of a 
solute. Considering the solvation of cylindric molecules of different widths and lengths, we 
study the penetration depth of the density field n(r) inside the solute molecule. 

The dependence of the solvation free energy AG cy \(L\R) upon its length L for two different 
radii R is shown in FigEk- 



IT 



< 




(a) 




FIG. 7: (a) Solvation free energy AG cy i(L|i?) of a cylinder as a function of its length L for different 
R in the range R = 0.9. ..1.5; (b) Slope of AG cy i(fi) as a function of R at L = 12 A. 

As one sees, from Fig|Hj3, the slope of the curves in Fig|7|becomes constant approximately 
for R> i? ~ 1.4f = 3.9 A. It means that for R > R the fluctuating density field penetrates 
inside the body of the solute to some fixed length only (to the depth of the "skin layer"). 
Hence, we may not account of all internal points of the solute located into its body beyond 
the skin layer Rq w 3.9 A. 

2. Secondary, we investigate the free energy of interaction AG sp h = AG sp hint(-D|-R) — 
AG sp hsoiv(-D — > oo\R) of two spherical solute molecules, each of radius R, upon the distance 
D between their centers. The free energy of interaction, or "the mean force potential" [32II33I] 
AG sph (D|_R), with the subtracted solvation energies of two solitary spheres, is computed for 
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the hard-core potential, U(D), where 

U{D) 



oo for D <2R 
for D > 2R 



i.e. we again take into account only the cavitation part of the free energy of interactions, 
neglecting the dispersion (attractive) part of the direct Van-der-Waals interactions between 
the impenetrable spheres of radius R. 

Even for such simplified situation we have found that the behavior of AG sp h(-D|-R) on 
D is qualitatively different for different values of R compared to £ [4q|]. Let us stress once 
more that AG sp h is the free energy of interaction between two spheres with the subtracted 
solvation energies of two solitary spheres AG sp h S oiv(-D — > oo\R). The observed dependencies 
are shown in FigJSJ 







| R=0.95 







D, [A] 
(b) 



D, [A] 

(o) 



FIG. 8: Free energy of interactions AG sp h of two spherical solute molecules as a function of the 
distance between their centers D for R = 0.4£ (a); R = 0.7£ (b) and R = 0.9£ (c). 

For all radii, R = 0.4£ A , 0.7£ A and 0.9£ A we reproduce the non-monotonic dependence 
AG sp h (-D) on D. This non-monotonicity is the manifestation of the oscillatory behavior of 
the water correlation function. Such behavior has been found earlier in direct numerical 
simulations consistent with the predictions of the "Scaled Particle Theory" (see j32j and 
the references therein), and on the basis of estimation of the solvent accessible surface area 
|33| . The account of the Van-der-Waals attractive tail in the potential U(D) can essentially 
change the behavior AG sp h(-D), however the non-monotonic behavior AG sp h(-D) will still 
hold. 

Comparing FigjHk and FiglSt, we see that the dependence AG sp h(-D) is sensitive (even 
qualitatively) to the size R of interacting spheres. The effective attraction between smaller 
spheres is weaker than that of larger spheres. This effect can be easily understood. For 
sufficiently large spheres, R > £, there is an penalty in the free energy of keeping two solute 
molecules close to each other, which is the "classical manifestation" of the Casimir-type 
effect. At the same time, the sufficiently small solute molecules (with R<£) are located in 
spontaneously created fluctuational cavities of size £ in water and the penalty for expelling 
the water from the volume occupied by two solute molecules close to each other is essentially 
smaller than that for large solutes. Moreover, one can see from FigJS^i that the mean force 
potential has the attractive (negative) part corresponding to the second water shell only. So, 
the two sufficiently small spherical solutes can form a bound state separated by one water 
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shell. The found behavior qualitatively coincides with the results of the work [32j, where the 
authors have studied the "potential of the mean force" acting between two spherical solute 
molecules. 



C. Electrostatics in fluctuating dipolar environment 



There are different ways to implement the electrostatic interactions into our self- 
consistent mean-field description. In many works it is assumed that the dielectric permit- 
tivity of the solvent, e, linearly depends on the solvent density, n(r), i.e., e(r) = 1 + an(r). 
In our approach we get rid of such supposition and follow the self-consistent scheme where 
the solvent is considered as a "gas" of dipoles. Our method ideologically is very close to 
the consideration of the Debye screening in electrically neutral plasma [34| . However, to the 
contrast with plasma, in our case the positive and negative ions are connected by holonomic 
constraints in pairs, forming the short-ranged dipoles. The screening in the system of ex- 
tended objects was a subject of many investigations. The description, most appropriate for 
our goals, has been developed in the appendix of the paper |3o| . 

Below we derive the basic set of equations describing simultaneously electrostatic and 
hydrophobic interactions in the fluctuating media. The detailed analysis of the obtained 
equations together with the numerical computation of the free energies of solvation and in- 
teractions of charged molecules, will be the subject of the forthcoming publication [T3|, while 
here we analyze briefly the limiting case, in which the proposed theory coincides with the 
electrostatics in the continuous media with some effective density-dependent permittivity 
£ (r). 

Let us introduce the quantities: 

Pmoi( r ' r ) ~~ the charge distribution of the solvent molecule with the orientation r and with 
the center of mass situated at the point r. Outside of the molecule Pmoi( r > r ) = 0> 

Psoiute( r ) ~~ the charge distribution of the solute; 

e fr( r ) - the effective mean electrostatic field. 

Sticking to the mean-field description developed in |35j], we can place each solvent 
molecule in a self-consistent field electrostatic field, <f)(r), averaged over the volume of this 
solvent molecule. So, we have: 

0(r,r) = J <P cS (r')pl ol (r'-r,T)dr' (38) 

The electrostatic potential </> e ff(r) is determined by the Poisson equation 

Atx 

A0 eff (r) = p« olutc (r) - /i(r) (39) 

^•solute 

where /i(r) is the media charge density 

Mr) = / / (e-^ (rV) - l) pL.(r - r', r)p(r')^V (40) 



21 



and p(r) = n(r) + co(r) is the solute density. Neglecting the short-range fluctuations, u;(r) 
we may write in the mean-field approximation 

Mr) = / / ('"^ " X ) PlJr ~ r '' T ^ r '^ d3r ' ^ 
The full Hamiltonian of the system reads now: 

if[n(r),&ff(r)] = ^ (^(Vn(r)) 2 + W») + i J dr dr'n(r) X ^(r, r')n(r') 

+ A / rfr(V0 cff (r)) 2 + / rfr^(r) (a^f) + -^PsliuteW + M')) 

J J V tsolute / 

(42) 

Let us stress that according to 36] we should minimize the functional action of hydrophobic 
interactions together with the one of electrostatic field, while the Poisson equation (|39j) ) 
is considered as a constraint and hence enters in the action with the functional Lagrange 
multiplier <^( r )- 

Minimizing the functional H[n(r), e ff( r )] with respect to the fields (p(r), </> e ff( r ), n(r), we 
come to the system of self-consistent Euler equations, fixing the corresponding equilibrium 
distributions. 

To simplify the current consideration, let us regard the case of point dipoles. The molec- 
ular density distribution p^ ol in this case is given by 

plJr) = ed(r - - e5 (r + ~ -e lV5(r) = -dW(r) (43) 

where d is the dipole moment of the molecule. Then 

0(r, d) = - J dr' cfr (r') dV r >5(r' - r) = dV0 cff (r) (44) 

this equation together with the linearized version of (|4f)|) leads to the following expression 
for the media charge density fi(r): 

Mr) = J J dr'^-pdV^ir') dV r 5(r-r') = ^ V (n(r)V0 cfr (r)) (45) 

Inserting the last relation into the Hamiltonian (|42j) and minimizing the corresponding 
expression, we get 

— = -oAn(r) + — + / tZr'Xta (r, r')n(r') - _£Vp(r)V&ff(r) = 
= d ll V (n(r)V V (r)) - A A</> cfr + Ay>(r) = (46) 

<50eff 3 V J 

— = A</> cfr + — pLutc(r) + -^V(n(r)V0 eff (r)J = 
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This system should be provided with some physically justified boundary conditions, for 
example 

► for |r| — > oo; e ff( r ) — * for |r| — ► oo; n(r) — > n for |r| — > oo (47) 

The numerical computation of the solvation and interaction free energies of charged ex- 
tended objects in the dipolar fluctuating media (in the water) based on the set of derived 



equations is in progress and will be the subject of our forthcoming publication [17 . 

It is easy to check that in the limit of the weak electrostatic field equations (|35Jl lead 
to the standard electrostatics in the continuous media with an effective dielectric density- 
dependent permittivity. Actually, when the field (p e g is small, we may assume 

if ~ eff (48) 

what means that 

VpV0eff~^ (49) 

We have omitted in (|3Sj) - (|35|l all numerical prefactors, supposing for simplicity that all 
quantities are dimensionless. 

Taking into account (pfS|)-eq:el9. we can rewrite the system (J41jJ) in the lowest order with 
respect to the electrostatic field: 

-aAn(r) + + f dr' X in ( r > r')n(r') = 

Jv * (50) 
v(£(n(r))V^(r)) + % s 9 olute (r) = 

where 

d 2 3 

e(n(r)) « 1 + -^-n(r) (51) 
3 

The closed set of equations (|5T)j) -()51 j) describes the electrostatics with the dielectric permit- 
tivity linearly proportional to the solvent density. 

Let us stress that our derivation implies the absence of the correlations between the 
orientational degrees of freedom of water molecules in the bulk and near the solute surface — 
only in this case we can pre-average the density /x(r) in (|4*T|) over the orientational degrees 
of freedom. We support our supposition by the conclusion made in |37l . 138^ on the basis of 
extensive numerical simulations: the rearrangement of the orientational degrees of freedom 
of water molecules near the solute surfaces plays the secondary role in the thermodynamics 
of solvation (see also for more discussions). 

To conclude with, let us repeat that in the current work the most attention is paid 
to the consideration of the pure hydrophobic (cavitation) effect. In particular, the main 
achievements are as follows. 

1. The method described in the paper provides the constructive basis for the consistent 
description of solvation effects of hard solute molecules of any shape. One sees that the con- 
tinuous two-length scale fluctuational approach with only two free parameters tuned to the 
Monte-Carlo results on solvation of hard spheres, manages to describe the known cavitation 
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contributions to the solvation free energy of alkanes without any additional adjustment of 
these parameters. 

2. In the framework of the same approach we have considered two auxiliary problems: 
a) the computation of the solvation free energy of solutes of cylindric shape, and b) the 
computation of the free energy of interactions of two spheres separated by some distance. 
The consideration of the problem (a) allows one to conclude that the depth of the "skin 
layer", i.e. the penetration length of the fluctuating field into the body of the solute is of 
order of 1.4£ A, where £ is the characteristic length scale of the theory (the location of the 
first maximum in the bulk correlation function of the solvent) while the investigation of the 
problem (b) permits us to make conjectures about the potential of the mean force acting 
between hard spheres in the fluctuating media. 

3. The developed approach can be extended to take into account the dispersion part of 
the solute-solvent interactions, as well as the electrostatic contributions of charged solute 
molecules. We describe the basic steps towards the implementation of both these contribu- 
tions into the current approach. 
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